
% save fig_trunc95.mat
%clear all
%load fig_Reiter.mat
%save fig_Reiter95.mat

clc;
clear;
%close all;
rho_u = 0.1;
load ../toDynare

omega_h = Ramsey.weights.omega_h;     %per history





name = ['fig_opti',num2str(rho_u),'.mat'];




load (name,'ut_eps','Yt_eps','Lt_eps','Ct_eps','Kt_eps','W_eps')

coef = 1;%ut_eps(1);

if eco.sigma==1
    util = @(c) log(c);
else
    util = @(c) (c^(1-eco.sigma)-1)/(1-eco.sigma);
end
trunc  = Ramsey.truncatedModel.truncatedAllocation;
x = trunc.c_h - (1/(eco.chi*(1/eco.phi + 1)))*(trunc.l_h).^(1+1/eco.phi);
S_h   = trunc.S_h;
W = sum(omega_h.*util(x));


%%% I construct the consumption equivalent for welfare
periods = length(ut_eps);
conso_eq1 = zeros(periods,1);

for i= 1:periods
    delta = 1;
    xtemp = trunc.c_h*delta - (1/(eco.chi*(1/eco.phi + 1)))*(trunc.l_h).^(1+1/eco.phi);
    Wtemp = sum(omega_h.*util(xtemp));
    while abs(Wtemp - (W+W_eps(i)))>1e-5
        delta = delta -0.1*(Wtemp - (W+W_eps(i)));
        xtemp = trunc.c_h*delta - (1/(eco.chi*(1/eco.phi + 1)))*(trunc.l_h).^(1+1/eco.phi);
        Wtemp = sum(omega_h.*util(xtemp));
    end
    conso_eq1(i) = (delta-1)*100; % percentage change
end

    

K1 = Kt_eps/coef;
Y1 = Yt_eps/coef;
L1 = Lt_eps/coef;
C1 = Ct_eps/coef;

ut1 = ut_eps;






%%

%name = ['fig_opti',num2str(rho_u2),'.mat'];
name = ['fig',num2str(rho_u),'FB.mat'];
load (name,'ut_eps','Yt_eps','Lt_eps','Ct_eps','Kt_eps','Ctot','Ltot','W_eps')

%coef = ut_eps(1);


%%% I construct the consumption equivalent for welfare
periods = length(ut_eps);
conso_eq2 = zeros(periods,1);
W = Ctot - (1/(eco.chi*(1/eco.phi + 1)))*Ltot.^(1+1/eco.phi);

for i= 1:periods
    delta = 1;
    xtemp = Ctot*delta - (1/(eco.chi*(1/eco.phi + 1)))*Ltot.^(1+1/eco.phi);
    Wtemp = sum(xtemp);
    while abs(Wtemp - (W+W_eps(i)))>1e-5
        delta = delta -0.1*(Wtemp - (W+W_eps(i)));
        xtemp = Ctot*delta - (1/(eco.chi*(1/eco.phi + 1)))*Ltot.^(1+1/eco.phi);
        Wtemp = sum(xtemp);
    end
    conso_eq2(i) = (delta-1)*100; % percentage change
end



K2 = Kt_eps/coef;
Y2 = Yt_eps/coef;
L2 = Lt_eps/coef;
C2 = Ct_eps/coef;

ut2 = ut_eps;





%%
LL = length(ut_eps);%length(B_eps);
LL = 200;
X=1:1:LL;

ligne = zeros(LL,1);

nx = 1;
ny = 5;

figure;
 
%subplot(nx,ny,1);
%plot(X,ut1(1:LL),'k-',X,ut2(1:LL),'b--','LineWidth',2);
%title('1. G','Interpreter','latex');
%ylim([-0.4 1.5]);


subplot(nx,ny,1);
plot(X,Y1(1:LL),'k-',X,Y2(1:LL),'b--','LineWidth',2);
title('2. $Y$','Interpreter','latex');
hold on
plot(X,ligne,'r--');
%ylim([-0.8 0.1]);
 

subplot(nx,ny,2);
plot(X,K1(1:LL),'k-',X,K2(1:LL),'b--','LineWidth',2);
title('3. $K$','Interpreter','latex');
legend('IM','FB','Location','southeast')
legend boxoff


subplot(nx,ny,3);
plot(X,L1(1:LL),'k-',X,L2(1:LL),'b--','LineWidth',2);
title('4. $L$','Interpreter','latex');
hold on
plot(X,ligne,'r--');
%ylim([-0.02 0.]);


subplot(nx,ny,4);
plot(X,C1(1:LL),'k-',X,C2(1:LL),'b--','LineWidth',2);
hold on
plot(X,ligne,'r--');
title('5. $C$','Interpreter','latex');
%ylim([-0.02 0.]);

subplot(nx,ny,5);
plot(X,conso_eq1(1:LL),'k-',X,conso_eq2(1:LL),'b--','LineWidth',2);
hold on
plot(X,ligne,'r--');
title('6. Welfare','Interpreter','latex');
%ylim([-0.02 0.]);





